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ABSTRACT 

We report contributions to cosmic infrared background (CIB) intensities originating from known 
galaxies and their faint companions at submillimeter wavelengths. Using the publicly-available Ul- 
traVISTA catalog, and maps at 250, 350, and 500 p.m from the Herschel Multi-tiered Extragalactic 
Survey (HerMES), we perform a novel measurement that exploits the fact that uncatalogued sources 
may bias stacked flux densities — particularly if the resolution of the image is poor — and inten¬ 
tionally smooth the images before stacking and summing intensities. By smoothing the maps we are 
capturing the contribution of faint (undetected in Ks ~ 23.4) sources that are physically associated, 
or correlated^ with the detected sources. We find that the cumulative CIB increases with increased 
smoothing, reaching 9.82 zb 0.78, 5.77 zb 0.43, and 2.32 zb 0.19 nWm“^sr“^ at 250, 350, and 500 p.m at 
300 arcsec full width at half-maximum. This corresponds to a fraction of the fiducial CIB of 0.94zb0.23, 

1.07 zb 0.31, and 0.97 zb 0.26 at 250, 350, and 500 [rm, where the uncertainties are dominated by those 
of the absolute CIB. We then propose, with a simple model combining parametric descriptions for 
stacked flux densities and stellar mass functions, that emission from galaxies with \og{M/M q) > 8.5 
can account for the most of the measured total intensities, and argue against contributions from ex¬ 
tended, diffuse emission. Finally, we discuss prospects for future survey instruments to improve the 
estimates of the absolute CIB levels, and observe any potentially remaining emission at z > 4. 

Subject headings: cosmology: observations, submillimeter: galaxies - infrared: galaxies - galaxies: 
evolution - large-scale structure of universe 
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1. INTRODUCTION 

Of all the light that has been emitted by stars, about 
half has been absorbed by interstellar dust and thermally 
re-radiated at far-infrared to submillimeter wavelengths, 
appearing as a diffuse, extragalactic, cosmic infrared 
background spanning 1-1000 Rm (CIB; iHauser fc DwekI 
12001^ [Dole et al.l 120061) . Statistically characterizing the 
sources responsible for this background is necessary to 
gain a full understanding of galaxy formation and cos¬ 
mology, and thus remains an ongoing pursuit. 
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The CIB was first detected in spectroscopy with 
the Far Infrared Absolute Spectroph otometer (FIRAS; 
iPuget et alJ 119961 : iMather et al.lIl999li . Observations of 
local starburst galaxies with the Infrared Astronomical 
Satellite 1/ii’AiS' nSoifer et all 1 19841 1 showed that galaxies 
could emit a surprisingly large part of their energy in the 
far-infrared, and ground-based measurements later con¬ 
firmed the existence of a rare population o f extremely 
lumin ous submillimeter galaxies (or SMGs; iBlain et alJ 
l2002ll . While these bright objects generated tremen¬ 
dous excitement, and in fact constitute a significant frac¬ 
tion of the total star-formation rate density at z > 3 
(e.g., iLe Floc’h et al.l [20051 iMuruhv et alJ [2?)TTl l. their 
low abundance only accounts for a small fraction of the 
total CIB. 

The arrival of the Herschel Space Observatory — 
whose instruments, PACS (70, 100, and 160 pm; 
Poglitsch et nilmoll and SPIRE (250, 350, and 500 pm; 


Griffin et al.ll2010[) 7 bracket the peak of the thermal spec¬ 
trum of dust emission — brought the promise of directly 
detecting less luminous and far more numerous dusty 
star-forming galaxies (DSFGs). However, source con¬ 
fusion resultin g from the relativel y large point-spread 
functions fe.g.. iNguven et ffi1l2010ll limited the number 
of galaxies that could be indi vidually resolved by P ACS 
at 10 0 and 160 pm to 75% (|Magnelli et al.l 1201311 and 
74% (|Berta et al.l 1201111 of the CIB, respectively, and 
by SPIRE a t 250 pm to 15% (iBethermin et al.l 120121 : 
lOliver et al.ll2012ll. Statistical methods including s tack- 
ing fe. g.. IBethermin et al.l|2012t IViero et al.ll2013all and 
P(D) (iGlenn et al.l l2010nBerta et al.l 1201111 performed 
better, resolving 70% at 250, 350, and 500 pm for the 
former, and 89% and 70% of the CIB at 100 and 250 pm 
for the latter, respectively. 

The origin of the rest remained unclear. IViero et al.l 
(|20I3all suggested that the missing flux could be tied 
up in faint sources — faint either because they are 
low mass, at high redshift, or extremely dusty. An¬ 
other possible source is diffuse emission from the dust 
that is known to be distributed in the halos of galaxies 
(|Menard et al.l[2010HHildebrandt et al.ll20I3ll. and could 

be heated by stripped s tars fe.g., iTal &: van DokkumI 
1201 It IZemcov et al.l 201411 . 

Here we present a technique to demonstrate that most, 
if not all, of the CIB can be accounted for by the com¬ 
bined emission from galaxies detected in current near- 
infrared surveys, and their faint companion objects, at 
z < 4. We show with a simple model that galaxies alone 
are the most plausible source of this signal, and argue 
that any remaining CIB likely originates from galaxies 
at still higher redshifts. 

2. DATA 

2.1. UltraVISTA Catalog 

We perform our analysis on catal ogs and im- 
ages located in the COSMOS field (IScoville et al.l 
[2^ . centered at 10^00“26",-f2°13'00". We use the 
Ks = 23.4 (AB)-selected , publicly-available^ catalog 
from iMuzzin et all (|2013bl . UltraVISTA), which consists 
of a 1.62 deg^ subset of the full COSMOS field, where 
both near-infrared and optical wavelengths are available 


(30 bands in all). The catalo g contains photometri c red- 
shifts computed with EAZY (|Bramm er et al.l l2008D. an d 
stellar masses computed with FAST ( Kriek et al.ll2009ll . 
Galaxies are split into star-forming or quiescent based 
on their positions in the rest-frame U — V vs . V — J 
color-color diagram 1 f/UJ: IWilliams et al.ll2009tl . 

2.2. Herschel/Her MBS Submillimeter Maps 

We use submillimeter maps observed with the 
Spectral and Pho tometric Imaging REceiver (SPIRE; 

I Griffin et al.l I20I0I1 at 250, 350, and 500 pm from the 
Herschel Multi-ti ered Extragalactic Survey (HerMES; 
lOliver et al.ll20I2ll . COSMOS is a level 5 field, consist¬ 
ing of 4 repeat observations, to an instrumental depth of 
15.9, 13.3, and 19.1 mJy (Scr), with confusion adding an 
additional noise term of 24.0, 27.5, and 30.5 mJy (5g) a t 
250, 350, and 500 pm, respectively ( Nguven et al.l 201^. 
Absolute calibration is detailed in iGriffin et al.l (l2013f) , 
with calibration uncertainties of 5%. Maps are made 
with 4 arcsec pixels at all wavelengths us ing the SMAP 
(|Levenson et al.l[201flHViero et al.ll2013bll pipeline. 

SPIRE maps are chosen specifically for this study 
because its wavelengths probe the rest-frame peak of 
therm al dust emission at z = 1-3 (iMadau fc DickinsonI 
I20I4I1 . and because large-scale features can be re con- 
structed with minimal filtering (|Pascale et al.l[2Mll) . 

3. METHOD 

We now present a novel method to estimate the 
extent to which known sources and their faint com¬ 
panions contribute to the CIB. We do this by ex¬ 
ploiting an inherent weakness of stacking: that stack¬ 
ing on images with poor angular resolution can re¬ 
sult in a boosted (or biased) average flux density 
arising fr om faint, uncatalogued, companion galax¬ 
ies ( e .g.. iSerieant et al.l 1200^ iFernandez-Conde et akl 

na 


20101: iKurczvnski &: Gawised 120101 : iHeinis et al.l 120131 : 


Viero et al.112013^ The trick is in recognizing that 


only correlated (i.e., clustered) sources will bias the 
stacked flux density (for an in-dept h discussion see 
iMarsden et al.l[^009t Iviero et al.ll2013all . so that summed 
intensities estimated with increasingly smoothed maps 
places limits on the full intensity in a given redshift range. 
Meanwhile, emission that is not correlated — say, emis¬ 
sion coming from sources at redshifts greater than those 
of our catalog objects — will not influence the primary 
measurement, except as a potential noise term. 

For this analysis we use the publicly-av ailable SIM- 
STACK code^, which is described in detail in IViero et al.l 
(|2013all . Briefly, synthetic images of the sky are con¬ 
structed from correlated subsets of catalogs (i.e., objects 
in the same redshift range) with the assumption that 
galaxies that are physically similar — in this case quies¬ 
cent or star-forming galaxies within a stellar mass bin — 
have comparable infrared luminosities and submillimeter 
flux densities. Synthetic images are then convolved with 
the PSF of the instrument, and are fit simultaneously to 
the actual sky map to retrieve the mean flux densities of 
the subsample. 

To induce a bias we smooth the maps by convolving 
them with Gaussians whose widths are the geometric dif¬ 
ferences of the nominal SPIRE and effective beams, i.e.. 


^ http://www.strw.leidenuniv.nl/galaxyevolutlon/ULTRAVISTA 


http://www.stanford.edu/-viero/downloads.html 
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Fig. 1. — Cumulative measured C IB vs. the size of t he effective image resolution (in arcsec, FWHM) at 250 (left panel), 350 (center 
panel), a nd 500 um (right pane l). The lFixsen et all d 19981) FIRAS values and 1 cr errors are shown as dashed line s with gray hatched regions, 
while the lLagache et ^ (119991) FIRAS values are shown as 3-dot-dashed lines, and the lBethermin et al.l (120121) model estimates are shown 
as dot-dashed lines. Colors represent the sum over all bins up to the given z. Grey dotted lines show the full set of null tests, and shaded 
regions the 1 u uncertainties, for the 2 < 4 case. The cumulative CIB vs. effective resolution increases more rapidly at higher redshift, 
where the catalog in increasingly incomplete. The flattening of the curve at the highest redshift suggests that any potential remaining 
intensity lies at higher redshifts. 


^smooth = \J o’gff — CTgpjj^g, where the nominal SPIRE 
resolutions are 17.5, 23.7, and 34.6 arcsec full width at 
half-maximum (FWHM) at 250, 350, and 500 nm, re¬ 
spectively, and the effective resolutions are 20, 30, 40, 
60, 80, no, 140, 180, 220, 260, and 300 arcsec FWHM. 
Above 300 arcsec, statistical uncertainties become con¬ 
sistent with zero (see Section a. All synthetic images 
and sky ma ps are mean - subtrac ted before stacking. 

Following IViero et al.l (l2013all we split the sample into 
8 x 8 bins of redshift (z = 0 to 4 with Az = 0.5) and 
stellar mass (5 star-forming and 3 quiescent), and stack 
the subsamples with SIMSTACK at each smoothing scale. 
Stacked flux densities are then color-corrected to account 
for the observed spectral shape of the sour ces through 
the pas sbands, with temperatures taken from iViero et all 
(|2013al Equation 18). Color corrections range between, 
1.006-0.984, 0.993-0.977, and 1.007-0.991 at 250, 350, 
and 500 pm, respectively. 

We use simulations to check that stacking and smooth¬ 
ing with a Gaussian, as opposed to the measured SPIRE 
beam or any other kernel, is a reasonable approximation, 
finding a bias of less than 4% for the largest smoothing 
kernel, which we include in the reported errors. We test 
that the method does not introduce unintended biases by 
performing 100 null tests for each set of synthetic images 
and maps. Null tests involve running the identical stack¬ 
ing pipeline with the same binning of sources, but after 
randomizing the positions of the sources in the catalogs. 
Because the map and images are mean-subtracted, we 
expect the stacked flux densities of the null tests to be 
consistent with zero. In Section 2] we show that this is 
the case. 

We note that this method has the advantage that miss¬ 
ing sources are not double-counted, meaning that the flux 
density from a single missing object will be distributed 
among the synthetic images, rather than appearing mul¬ 
tiple times (as would be the case in thumbnail stacking). 
Also note that stacked flux densities are intentionally not 


corrected for completeness because it is precisely the flux 
densities of the missing (i.e., incomplete, but similarly 
applies to misclassified AGN or DSFGs) sources that we 
are attempting to measure by degrading the effective res¬ 
olution of the map. As a result, this technique is not 
limited to GIB studies, but is applicable to any study 
where estimates of the level of faint, correlated emission 
are in question. 

4. RESULTS 

We report total intensities of our stacking measure¬ 
ment as the cumulative sum over redshift (from z = 0 
to the redshifts labeled with different color lines) vs. 
the effective resolution of the image, in Figure [Tf. At 
lower redshifts (z < 1) the fractional CIB that is mea¬ 
sured increases weakly with increasing effective beam 
size, which is expected given that the completeness of the 
catalog at these redshifts is near unity for stellar masses 
log(M/M 0 ) > 9. Conversely, at higher redshifts the 
fractional CIB that is measured increases rapidly with 
increased smoothing. The maximum CIB we resolve is 
9.82 ±0.78, 5.77 ±0.43, 2.32 ± 0.19 nWm-^sr-i at 250, 
350, and 500 pm. We note slightly different behaviors 
between the three bands; particularly at 500 pm, which 
appears to have not converged, and may be indicative of 
higher redshift contributions to the CIB. 

Estimates of the exact measured fraction of the abso¬ 
lute CIB are limited by the uncertainties in the reported 
ab solute values of the CIB derived from FIRAS spectra 
by iFixsen et al.l ((199811 and iLaeache et al.l (1199911 . which 
range between 22 % and 30%. The IFixsen et alJ (|1998[ 1 
levels (10.4 ± 2.3, 5.4 ± 1.6, 2.6 ± 0.6; hereafter chosen to 
represent the fidu cial values) are shown as dashed lines in 
FigurelU and the lLagache et al.l (Il999ll levels (11.8±2.9, 
6.4 ± 1 . 6 , 2.7 ±0.7) are shown as 3-dot-dashed lines. Ad- 

^ Tabulated values for the intensities in Figure 1 can be found 
at https://web.Stanford.edu/~viero/downloads.html 
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Fig. 2.— Left three panels: Stacked flux densities (mea sured with nominal PSFs) in divisions of stellar mass vs. redshift for star-forming 
galaxies. Power-law model fits (described in Section |5.Il l are shown as dotte d lines. Right panel: Parametric mo del for the s tellar mass 
functions of star-forming galaxies, which combines the parameterization of the lTomczak et al.l 1120141 . circl es! data bv |Leia et all II201SI . solid 
lines) at z < 2.5, and at modification at z > 2.5 (dashed lines) which interpolates to measurements from IMuzzin et aTTil^Mfl^ exes). 


ditionally, iBethermin et al.l (I2012L dot-dashed lines) pro¬ 
vide estimates of the total CIB by extrapolating their 
measured counts, finding they agreed with FIRAS with 
similarly large uncertainties . In to tal we find that the 
fraction of the iFixsen et al.l (I1998D fiducial CIB we re¬ 
solve is 0.94 ± 0.23, 1.07 ± 0.31, and 0.97 ± 0.26 at 250, 
350, and 500 pm, respectively. 

The full set of 100 null tests at each effective beam size 
are shown as gray dotted lines in Figure [TJ with the 1 a 
limits represented by the shaded grey region. For clarity 
only the tests for z < 4 are shown. As expected, the 
results of the null tests fluctuate around zero, with the 
magnitude of that fluctuation increasing with increasing 
effective beam size. Additional systematic uncertainties 
include calibration and beam area uncertai nties, and cos- 
mic variance, which is estimated following iMoster et al.l 

(Imh) . 

5. DISCUSSION 

The CIB can be divided into three components: (i) the 
contribution from the sources in the parent Ks ^ 23.4 

catalog; (ii) the contribution from correlated companion 
sources that for some reason do not make the catalog cut 
but are recovered with our smoothing-stacking method; 
and (iii) uncorrelated emission, including and likely dom¬ 
inated by sources lying at z > 4. 

It is interesting to consider the nature of (ii), the unde¬ 
tected sources that are captured by our smoothing proce¬ 
dure: are they very low-mass galaxies; massive but dusty 
galaxies that evade detection in current deep AT-selected 
catalogs; or some other unknown component? 

We now propose — through a simple model that com¬ 
bines a parametric description of the stellar mass func¬ 
tion with simple fits to nominal stacked flux densities — 
that the galaxies that make up the stellar mass function 
are alone able to describe our measurement, and that by 
extension the missing CIB can reasonably be attributed 
to galaxies in the low-mass end of the stellar mass func¬ 
tion. 

5.1. A Model of the CIB 

The fir s t com ponent of the model adopts the 
iLeia et ^ (|2015l Equation 1) parameterization of the 



Fig. 3.— Measurements (open circles) and models (solid lines) 
for cumulative CIB intensities vs. redshift. Measurements are the 
cumulative sums for images stacked at the highest effective b eam 
size of 300 arcsec FWHM. The models, described in Section 15.11 
combine parametric descriptions for the stacked flux densities (at 
the native resolution of the images) and t he stellar mass fu n ctions 
(see Figure [^. Extrapolated counts from l^thermin et al.l 1120121) 
are shown as gray squares. Also shown are cumulative CIB intensi¬ 
ties vs. redshift on stacks made with the native (i.e., non-smoothed) 
SPIRE images (diamonds), and the models after they have been 
modified to reflect the incompleteness of the actual catalog (dashed 
lines). 

iTomczak et al.l (l2014l) stellar mass function. The right- 
mos t panel of Figure [2] illustrates the performance of 
the iLeia et al. (I2015D model (solid lines) against the 
ITomczak et al (120141) data (circles). We find that it di¬ 
verges from the measurements at higher redshifts and so 
we add a modification at 2 > 2.5 (dotted lines) by simply 
interpolating betwee n the model and the m easured stel¬ 
lar mass functions of iMuzzin et al.l (|2013aD to 2 = 4. We 
check that the exact value of faint-end slope at high-2 
negligibly affects the integral, and set it to -1.6. 

Similarly, the second component of the model con¬ 
sists of power-law fits to the stacked flux densities vs. 
lookback time for each stellar mass bin, with the added 
condition that the mass dependence of the slopes and 
offsets themselves follow smooth functions. The three 
left panels of Figure [2] show the stacked flux densities 
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and best-fit power laws as open circles with error bars, 
and dotted lines, respectively. Finally, the model is in¬ 
tegrated over the redshift range z = 0.1-4, and stellar 
mass range log(M/M0) = 8-14; although the contribu¬ 
tion from galaxies below log(M/MQ) < 8.5 is found to 
be negligible. 

Figure |3] compares the resulting model (solid lines) 
with the full CIB measurements (open circles) at all three 
wavelengths. Note that model is not a fit, and yet de¬ 
scribes the measurement remarkably well, demonstrat¬ 
ing that the faint-end of the mass function can plausibly 
exp lain the recovered CIB. Also in good agreement are 
the iBethermin et al.l (j2012l gray squares) estimates de¬ 
rived through extrapolation of their number counts. We 
note s ome tension with the iPlanck Collaboration et ^ 
(j2014l) intensity at 350 p.m (7.7 ± 0.2nWm“^), which is 
an output of their halo model fit to CIB power spectra; 
although their 550 pm estimate from the same model is 
in relative agreement (2.3 ± 0.1 nWm“^). 

Also shown (as diamonds) are measurements made 
when stacking on the native SPIRE images (i.e., no 
smoothing). They are again compared to the model 
(dashed lines), but this time the model is modified so 
that the low-mass limit of the integral over the mass 
function effectively begins at higher masses with increas¬ 
ing redshift, reflecting the completeness behavior of the 
catalog. The difference between the solid and dashed 
lines can be interpreted as a measure of how much of the 
CIB originates from the parent sample, and how much is 
from companion sources not detected at this Ks limit. 
If, for arguments sake, the catalog were 100% complete 
at all stellar masses and all redshifts, then the dashed 
and solid lines would overlap. 

Arguments for additional, diffuse, sources of CIB; 
in particular dust in the extended halos of galaxies 
(iMfoard et al.l[2M0) . are thus disfavored. 

5.2. A CIB beyond z of 4? 

While our measurements are consistent with the fidu¬ 
cial levels of the total CIB, the existing uncertainties on 
the absolute level are so large that it is difficult to con¬ 
vincingly estimate how much CIB is still unresolved. Any 
missing CIB is more likely to occur at longer wavelengths, 
which are more sensitive to higher redshift s (because of 
the negative AT-correction JBlain et al.ll2003D . from where 
we would expect uncorrelated emission to originate. 

Indeed, several exceptional ULIRC-like galaxies have 
been identifie d at z > 4 , albe i t with low abun¬ 

dances f^e.g.. iRiechers et'all I20R^ IVieira et al.l 120131 
iDowell et ai.ll20l4 ISwinbank et al.ll2014[) . However, ex¬ 
trapolations of the contribution of ULIRGs to the star- 
formation rate den sity at z > 4 fe.g.. lMurDhv et aT1l2011l : 
IViero et al.l l2013all points to them dominating the far- 
infrared emission at this epoch. Determining the relative 
levels that they and less Inminous galaxies contribnte is 
a key qnestion going forward. 

To resolve these high-z questions, more data are re¬ 
quired; particularly: i) an npdate of the absolute CIB; 
ii) deeper catalogs with stellar masses and redshifts to 
redshifts greater than 4; iii) submillimeter surveys (350 


to 1000 pm) with large angular-scale fidelity and smaller 
PSFs in the regions of those deep catalogs — particularly 
at longer wavelengths which are more sensitive to galax¬ 
ies at higher redshifts. The former can only be achieved 
with space-based missions to measure the DC level above 
atmospheric foregrounds. The second requirement is 
steadily growing from multiple current or upcoming ef¬ 
forts (e.g., SDSS/BOSS, CANDLES, DES, LSST), al¬ 
though the photometric redshif ts for very dusty ga laxies 
may remain quite uncertain Isee lSpitler et al1l2014f) . The 
last point will require a ground-bas ed submillimeter ob - 
servatory similar in scope to CCAT (|Sebring et al.ll200^ . 

6. CONCLUSION 

We find that most of the CIB at 250, 350, and 500 pm 
can be accounted for by galaxies detected in current near- 
infrared surveys of moderate depth {Kg ~ 23.4), and 
galaxies correlated with them. We report total intensities 
of 9.82 ± 0.78, 5.77 ± 0.43, 2.32 ± 0.19nWm-2sr-i at 
250, 3 50, and 500 pm, wh ich corresponds to a fraction 
of the iFixsen et al.l (Il998f) absolute CIB of 0.94 ± 0.23, 
1.07 ±0.31, and 0.97 ±0.26. 

We find that a simple model combining parametric de¬ 
scriptions for the stellar mass function and stacked flux 
densities for log(M/M0) > 8.5 and z < 4 is able to con¬ 
vincingly reproduce the measurements, which supports 
the argument that the sources in the faint-end of the 
mass function make up the previously missing CIB. We 
note that emission from objects that are not in the cat¬ 
alog, and is either uncorrelated with catalog objects or 
exists at scales greater than 300 arcsec, would be missed. 
However unlikely, this cannot be ruled out without a bet¬ 
ter absolute measurement to compare against. 

Finally, we propose that any remaining CIB likely orig¬ 
inates from galaxies at z > 4, and if so should be de¬ 
tectable at submillimeter and millimeter wavelengths. 
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